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Is it appropriate to model turbidity currents with the 
three-equation model? 

Peng Thomas Pahtz,^’^ and Zhiguo He^’^ 

Abstract. 

The three-equation model (TEM) was developed in the 1980s to model turbidity cur¬ 
rents (TCs) and has been widely used ever since. However, its physical justification was 
questioned because self-accelerating TCs simulated with the steady TEM seemed to vi¬ 
olate the turbulent kinetic energy balance. This violation was considered as a result of 
very strong sediment erosion that consumes more turbulent kinetic energy than is pro¬ 
duced. To confine bed erosion and thus remedy this issue, the four-equation model (FEM) 
was introduced by assuming a proportionality between the bed shear stress and the tur¬ 
bulent kinetic energy. Here we analytically proof that self-accelerating TCs simulated with 
the original steady TEM actually never violate the turbulent kinetic energy balance, pro¬ 
vided that the bed drag coefficient is not unrealistically low. We find that stronger bed 
erosion, surprisingly, leads to more production of turbulent kinetic energy due to con¬ 
version of potential energy of eroded material into kinetic energy of the current. Fur¬ 
thermore, we analytically show that, for asymptotically supercritical flow conditions, the 
original steady TEM always produces self-accelerating TCs if the upstream boundary 
conditions (“ignition” values) are chosen appropriately, while it never does so for asymp¬ 
totically subcritical flow conditions. We numerically show that our novel method to ob¬ 
tain the ignition values even works for Richardson numbers very near to unity. Our study 
also includes a comparison of the TEM and FEM closures for the bed shear stress to 
simulation data of a coupled Large Eddy and Discrete Element Model of sediment trans¬ 
port in water, which suggests that the TEM closure might be more realistic than the 
FEM closure. 


1. Introduction 

Turbidity currents (TCs) are sediment-water mixtures 
rapidly-moving downslope through clear water. They sig¬ 
nificantly contribute to the evolutions of a large variety of 
sedimentary structures and morphological features in river 
reservoirs, lakes, estuaries, and deep oceans [Islam et al., 
2008; Meiburg and Kneller, 2010; Liu et al, 2012; Konsoer 
et al, 2013] and are thus of high interest for many fields of 
Earth Science. However, it has turned out quite difficult to 
carry out controlled in-situ measurements of TCs, explain¬ 
ing why only very few have been reported [Xu et al., 2004; 
Cossu and Wells, 2010; Pyles et al., 2013; Sumner et al., 
2013]. This highlights the importance of laboratory mea¬ 
surements and numerical modeling of TCs for developing a 
better understanding of their nature. 

There are two groups of numerical models: depth¬ 
resolving models [Strauss and Glinsky, 2012; Yeh et al., 
2013] and layer-averaged models, which include the three- 
equation model (TEM) and its variants [Fukushima et al., 
1985; Parker et al., 1986; Zeng and Lowe, 1997; Choi, 1998; 
Imran et al., 1998; Bradford and Katopodes, 1999; Kostic 
and Parker, 2006, 2007; de Luna et al., 2009; Toniolo, 2009; 
Hu and Cao, 2009; Kostic et al., 2010; Eke et al., 2011; Hu 
et al., 2012; Lai and Wu, 2013; Kostic, 2014; Elfimov and 
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Khakzad, 2014] and the four-equation-model (FEM) and its 
variants [Fukushima et al., 1985; Parker et al, 1986; Sala- 
heldin et al., 2000; Pratson et al, 2001; Das et at, 2004; 
Fildani et al., 2006; Kostic and Parker, 2006; Yi and Im¬ 
ran, 2006; Eke et al., 2011; Kostic, 2011; Tracer et al., 2012]. 
The FEM differs from the TEM in the way in which the 
bed shear velocity (u*) is computed [Fukushima et al, 1985; 
Parker et al., 1986]: while the TEM computes u* from the 
drag exerted on the bed, roughly approximated by 

TEM : ul = CdU^, (1) 

where Co > 0 is the bed drag coefficient and U the layer- 
averaged velocity of the sediment-water mixture, the FEM 
computes u* from the assumption that the bed shear stress 
is proportional to the layer-averaged turbulent kinetic en¬ 
ergy (k), 

EEM : ul = ak, (2) 

where a > 0 is the dimensionless proportionality constant. 
The inclusion of k in the FEM makes it necessary to also 
include the turbulent kinetic energy balance. This explains 
the different names of the models, which suggest that the 
EEM contains one governing equation more than the TEM. 
However, in our opinion, these names are slightly misleading 
since the same turbulent kinetic energy balance can also be 
used in TEM to compute k [Fukushima et al, 1985; Parker 
et al, 1986]. The actual difference is that k influences the 
evolution of TCs in the FEM, while it does not do so in the 
TEM. 

Why were two kinds of layer-averaged models, the TEM 
and FEM, developed for TCs? The answer is that the 
steady TEM simulations by Fukushima et al. (1985] (F85) 
and Parker et al. (1986] (P86) failed to reproduce physi¬ 
cally realistic self-accelerating TCs, which occur when the 
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bed slope (S) is sufficiently large to ensure that U and the 
sediment transport rate (■0) increase downstream without 
limit ever after a sufficiently large, finite distance down¬ 
stream {U,ip) - > ( 00 , 00 ), where x is the stream- 

wise coordinate). In fact, the authors found that the di¬ 
mensionless net production rate of the turbulent kinetic en¬ 
ergy {AE = {h/U^)dk/dx, where h is the depth of the TC), 
composed of production through conversion from potential 
energy and dissipation through erosion and suspension of 
bed sediment, becomes negative when s —> 00 for self- 
accelerating TCs. It follows that the TCs should die, which 
is inconsistent with its self-accelerating property [Fukushima 
et ai, 1985; Parker et ai, 1986]. The authors linked this os¬ 
tensible failure of the steady TEM to a possible overestima¬ 
tion of the bed sediment erosion rate, caused by a possible 
overestimation of u* in Eq. (1), and remedied this issue 
by assuming Eq. (2), which limits m* through a simplified 
first-order relation with k. 

In this paper, we report simulations of self-accelerating 
TCs using the steady TEM and parameter values and em¬ 
pirical relations exactly as reported by F85 and P86. We find 
that the simulations by P86 do not produce self-accelerating, 
but instead decelerating TCs when the upstream boundary 
conditions (“ignition values”) specified by P86 are used. The 
authors obtained these ignition values by setting h{0) = 2m, 
and following the procedure by Parker [1982]. However, if 
instead h(0) = Im is used, the simulations by P86 do result 
in physically realistic self-accelerating TCs, in contrast to 
the claim made in this study that such TCs do not occur 
for the specified parameter range. Moreover, we also find 
that simulations by F85 do result in physically realistic self- 
accelerating TCs even for the same ignition values, in con¬ 
trast to the claim made in this study. In fact, we find that 
the downstream profile of dk/dx is actually exactly opposite 
to the one described in F85, and thus aI: > 0 (where the 
arrow denotes hereafter the limit x —>■ 00 , ^ = iima;_>oo •)> 
consistent with its self-accelerating property, indicating a 
possible error in the computations by F85. We sup port 
this claim with an analytical proof showing that AE > 0 
for self-accelerating TCs simulated with the steady TEM 
if Cd or alternatively S are not unrealistically low. This 
proof contains the derivation of the asymptotic behaviors 
(x —>■ 00 ) of quantities characterizing a self-accelerating TC, 
such as the Richardson number (Ri), which we show to be 
asymptotically constant. Erom an analytica l st ability anal¬ 
ysis of the steady TEM, we then find that Ri < 1 (super¬ 
critical flow, i.e., the densimetric Fronde number Fr > 1 
since Ri = 1/Fr^) is a necessary and sufficient condition for 
the existence of self-accelerating TCs. On basis of this re¬ 
sult, we provide a novel method to find ignition values which 
always result in self-accelerating TCs for simulations using 
the steady TEM. Another interesting finding of our study is 
that, surprisingly, a larger bed sediment erosion rate, Eg, re¬ 
sults in larger values of AE, even though the erosion of bed 
sediment dissipates turbulent kinetic energy. This is because 
eroded bed sediment increases the sediment mass and thus 
potential energy of the TC, which is then converted into tur¬ 
bulent kinetic energy downslope. Our study also includes a 
comparison of the TEM and FEM closures for the bed shear 
stress (Eqs. (1) and (2)) to simulation data of a coupled 
Large Eddy and Discrete Element Model of sediment trans¬ 
port in water [Furbish and Schmeeckle, 2013; Schmeeckle, 
2014], which suggests that the TEM closure might be more 
realistic than the FEM closure. 

In the following, we first briefly review the conservation 
equations governing the steady TEM and FEM as repo rted 
by F85 and P86 in Section 2. Then we proof that AS > 0 
if Cm and S are not unrealistically low for self-accelerating 
TCs simulated with the steady TEM in Section 3. This 


section also contains the proof that < 1 is a necessary 
and sufficient condition for the existence of self-accelerating 
TCs. Afterwards in Section 4, we present our simulations 
using the steady TEM and parameter values exactly as re¬ 
ported in F85 and P86 and show that these simulations are 
consistent with our analytical proof. There we also present 
a method to find ignition values which always result in self- 
accelerating TCs for simulations using the steady TEM and 
show that stronger erosion leads to larger positive values of 
AE. The latter is then explained in Section 5, which also 
includes the comparison of the TEM and FEM closures for 
the bed shear stress to the aforementioned numerical data, 
and conclude in Section 6. 


2. Conservation equations 


Although unsteady models might be more realistic, we 
here consider steady TCs (d/dt = 0) in order to be consis¬ 
tent with the original studies by F85 and P86. For this case, 
the mass conservations of the sediment-water mixture (Eq. 

(3) ) and the sediment carried by the current (Eq. (5)), the 
momentum conservation of the sediment-water mixture (Eq. 

(4) ), and the turbulent kinetic energy conservation of the 
sediment-water mixture (Eq. (6)) are written as [Fukushima 
et ai, 1985; Parker et ai, 1986] 


dh 

dx 

U~d^ 
h dxj) 
dx 
h dk 
IPdx 


—RiS + em(2 — O.bRi) + ^ + ^RiR^ 
1 — Ri 

RiS — euj(l -I- 0.5Ri) — ^ — ^RiR^ 
1-Ri ’ 



-e^(l — Ri) + 



ke.g, 

IP 


r..Vs 


1 


RiR^ ) 0(fc) = AE, 


(3) 

(4) 

(5) 

( 6 ) 


where 0 denotes the Heaviside function, Ri = gip/U^ the 
Richardson number with g = p(ps — p'w)lpw > 0 the sub¬ 
merged value of the gravity constant (g) and pa (pw) the 
density of sediment (water), Us > 0 is the sediment settling 
velocity. 


ipe = EghU/ro 


(7) 


is the equilibrium sediment transport rate at which sedi¬ 
ment exchange between the current and the bed vanishes 
(Rp = 0), while Bu, > 0 (water entrainment rate). To > 1 
(ratio between bed and average sediment concentration). 
Eg > 0 (bed sediment erosion rate), and /3 = heojt?^'^ > 0 
(co > 0 is the average rate of viscous dissipation of turbu¬ 
lent kinetic energy) are bounded coefficients. We note that 
Eq. (6) has been slightly modified from the version reported 
by F85 and P86, namely it has been multiplied by 0(fc). 
While this modihcation has no relevance for practical appli¬ 
cations because 0(fc) = 1 if fc > 0, we incorporated it here 
for the mathematical proof in Section 3 since it ensures that 
AE does not become a complex number due to 0(fc) = 0 if 
fc < 0. Indeed, without 0(fc), fc could become negative and 
thus AE complex due to the term fc^^^. Eqs. (3-5) consti¬ 
tute the steady TEM together with the closure Eq. (1). In 
contrast, the steady FEM is constituted by Eqs. (3-6) since 
the closure Eq. (2) incorporates fc, which must be computed 
by Eq. (6). However, even though Eq. (6) does not influence 
the values of h, U, and -ip in the steady TEM, it is still used 
to compute AE if required [Fukushima et al., 1985; Parker 
et al, 1986]. Moreover, it is important to point out the 
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fact that the computed layer-averaged volumetric sediment 
concentration, defined by 



( 8 ) 


must be smaller than unity, which is automatically ensured 
by the steady TEM for most practical applications (mainly 
due to Eq. (5), which makes ip decrease strongly when C 
becomes large). In fact, as we show in Section 3, Eq. (8) is 
always fulfilled for self-accelerating TCs in the limit x ^ oo. 

F85 and P86 used the following empirical relationships to 
compute the coefficients e™, Xo, Ea in the steady TEM and 
FEM and P in the FEM, 


To 

Ea 

P 


0.00153/(0.0204 -bffi), 
l + 31.5(u./u,)"^-'‘®, 

[ 0.3 

\ 3 X 10-i2Zi°(l - h/Z) 
0 


^ > 13.2 
5 < ^ < 13.2 
^ < 5 


(9) 

( 10 ) 

( 11 ) 


0.5e™(l -Ri- 2Co/a) -b Co 
{Co/aY-^ 


(only FEM), (12) 


where Z = y/Rep{ut/va) with Rep = ^/gD^ji/ the particle 
Reynolds number, Da the mean sediment particle diameter, 
and v the kinematic viscosity of clear water. From simula¬ 
tions with the stea dy TEM using Eqs. (9-11) and S = 0.08, 
F85 obtained < 0 for their simulated self-accelerating 
TCs and thus concluded that the TEM produces physically 
unrealistic results. In the following section, we anal ytically 
show that the steady TEM can only result in < 0 if 
Cd and S are both smaller than certain threshold values 
defined later. In particular, if Eq. (9) is used to compute 
e-w, Cd must be smaller than 0.00097 (which is much smaller 
than the authors’ Cd = 0.004) and S be smaller than 0.0073 
(which is much smaller than the authors’ S = 0.08). 


3. Analytical proof 


the property U = ip = Our proof does not require 
particular empirical expressions or values for the empirical 
parameters Cd, e™. To, Ea, and P, such as Eqs. (9-12). 
It only requires that these parameters are bounded as well 
as lim( 7 _>oo Ea > Q (bed material must be eroded if U is 

sufficiently large), e^iRi < oo) > 0, Cu, — ^~^°°> 0 (water 
is entrained if and only if flow turbulence is present), and 
de-w/dRi < 0 (water entrainment increases with turbulence). 
All properties needed for the proof are summarized in Ta¬ 
ble 1. 

Moreover, in our proof we often formally operate with 
quantities in the limit a; —>■ oo, which requires that this limit 
exists for these quantities (note that a limit also exists, if it is 
infinite). However, functions with spatial periodicity might 
not fulfill this requirement (e.g., sin® does not exist). Since 
Eqs. (3-6) do not explicitly contain periodic functions, it is 
safe to presume that such limits always exist for our case. 
Finally, we will often use the following rewritten versions of 


Eq. (5), 



Rtp = 

EaVah/lp — XoVa/U, 

(13) 

Ft-ip — 

WaVah/lp, 

(14) 

dip/dx = 

EaVa — XoVaC, 

(15) 

dip/dx = 

~^aVa — rtva~d. 

(16) 


where we used Eqs. (7) and (8), and U = oo, and that 
Ea > 0 (from lim(7_>oo Ea > 0), Ea < oo, 0 < Va < oo, and 
0 < rt < 00 . Eqs. (14) and (16) are the limits ® —>■ oo of 
Eqs. (13) and (15), respectively. 

Our pyoofps separated into three parts. First, we cal¬ 
culate Ri, and the asymptotic behaviors of U, ip, and 
h in Section 3.1. Using th e r esults of Section 3.1, we then 
proof in Section 3.2 that Ri < 1 is a necessary and suf¬ 
ficient condition for the existence of self-accelerating TCs, 
and in Section 3.3 th at P d > Comin or S > Rmin are suf¬ 
ficient conditions for AE > 0. Afterwards we discuss why 
these condition are virtually always fulfilled for physically 
relevant cases in Section 3.4. 


In this section, we proof that Ai? > 0 for self-accelerating 
TCs simulated with the steady TEM if Co > Comin or 
S > Smin, where Comin and Rmin are certain values of Co 
and respectively, which we define later. We also proof 
that Ri < 1 is a necessary and sufficient condition for the 
existence of self-accelerating TCs. For the proof, we make 
use of the definition of self-accelerating TCs, which includes 


Table 1. Summary of properties needed for the analytical 
proof 


To show: l} = ij) = oo => A.^ > 0 if Cd > C'omin or S' > Smin 


Property 

Explanation 

Eqs. (1) and (3-6) 

definition of the TEM 

—>■ A -4 

h>0, f/ = '0=OO 

definition of self-accelerating current 

0 < 9 = 9 < oo 

constant, positive parameter 

0 < S = S < oo 

constant, positive parameter 

0 < Cd = Cd < oo 

constant, positive parameter 

0 < Va = vt < oo 

constant, positive parameter 

1 < rt < OO 

To > 1 is a bounded parameter 

0 < /3 < oo 

/3 > 0 is a bounded parameter 

0 < ef/ < cxD 

e^t; > 0 is a bounded parameter 

0 < Eg < oo 

Eg > 0 is a bounded parameter 

limj7_>oo Ea > 0 

erosion occurs for infinite current speed 

ew{Ri < oo) > 0 

water is entrained by turbulence 

— 0 

no entrainment without turbulence 

dewjdRi < 0 

entrainment increases with turbulence 


3.1. Asymptotic solutions 

This part of the proof follows several logically ordered 
steps to show 0 < Ri < 00 and R^ < oo. (Note that 
Ri > 0 is not trivial, even though Ri > 0 is fulfilled, since 
in the limit this quantity could approach the boundaries of 

its restricted domain.) Afterwards we use these results to 
-^ ^ 

calculate R^, Ri and the asymptotic behaviors of U, ip, and 
h. Note that all results of this sectio n are for mally also valid 
for the FEM if one replac es Cp^ by ak/U^ (from Eqs. (1) 
and (2)) and assumes 0 < ak/U^ < oo. 

3.1.1. Showing Ri > 0 

Ri > 0 can be shown by presuming Ri = 0 and arriving 
at a contradiction. Inserting this presumption in the limit 
® —>■ oo of Eq. (4) using Eq. (1) yields 

< 0, (17) 

U ax 2 

where we used eft > 0, Cd > 0, and Ril^ > 0 (from 
if = ^ = oo). Eq. (17) is a contradiction to 1/ = oo, which 
requires dU/dx > 0 and thus A ^ > q. Hence, Ri > 0. 

3.1.2. Showing ^ < oo 

While it is physically clear that C cannot be larger than 
unity, we do not need to assume this beforehand for the ana¬ 
lytical proof. Instead, this property is strictly obtained from 
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the model equations and properties of self-accelerating TCs. 
Here we first show that o < oo, while we later even obtain 
C? = 0 when calculating the asymptotic profiles. C? < oo 
can be shown by presuming ~^ = oo and arriving at a con¬ 
tradiction. From C? = 00 follows that 


d’^jdx = '^sVs — 7 ^Vb~^ = — 


oo, 


(18) 


where we used Q < Es < oo, 0 < Vs < oo, > 0, and 
Eq. (16). Eq. (18) is a contradiction to i/) = oo. Hence, 
~d <oo. 

3.1.3. Showing h — oo 

h = oo can be shown by presuming h = he, where he 
is a finite length, and arriving at a contradiction. Inserting 
this presumption in the limit a; —>■ oo of Eq. (8) yields 


where we used Ri < oo and sgn denotes the slgnum func- 
tion. Eqs. (26) and (27) mean that either dh/dx or E ^ 

will become —oo depending on the limit of th e sign of 1 — Ri. 
However, a n^ative value of dh/dx or is a contradic¬ 

tion tolj — h = 00 . Hence, 

3.1.6. Asymptotic behaviors of U, l/>, and h 

The asymptotic behavi ors of U, tp, and h can be calcu¬ 
lated through computing R^ and Ri using I’Hospital’s rule 
[Chatterjee, 2012). First, if R^ > 0, using R^ < oo, we can 

use the same arguments which we used prior to Eq. (25) 
and calculate R^ by 


^ = -^RiU^ = 00 , 

he heQ 


(19) 


where we used ~q < oo, Ri > 0, and ij = oo. This is a 

—>■ 

contradiction to C < oo. Hence, h = oo. 

3.1.4. Showing Ri < oo and > 0 

Ri < oo can be shown by presuming Ri = oo and arriving 
at a contradiction. Inserting this presumption in the limit 
a; —>■ 00 of Eqs. (3) and (4), respectively, yields 


->• 


dh/dx = S —-R^, 

yW 

U dx 


= —S + = —dh/dx - 


w 


h 


U dx 


+ dh/dx = 0 


( 20 ) 

( 21 ) 

( 22 ) 


where we used S < oo, Cd < oo, Eq. (1), an d e^, g. 

Due to h = 00 and l) = oo, dh/dx > 0 and E ^ > g must 
be fulfilled. It then follows from Eq. (22) that 


~~dxf 


dh/dx = 77 -;— = 0 

U dx 


(23) 


and thus R^ = 25 due to Eq. (20). Since 0 < 5 < oo, this 
allows us to calculate 


dip/dx = R.,pp}/h = R,/, ip/h = R,/, - y = EsVs, 

h/ip 


(24) 


where we used Eqs. (5) and (14). Using this and I’Hospital’s 
rule [Chatterjee, 2012], we then also obtain from Eq. (14) 


0 < 25 = R,/, = EaVsh/p) = EaV, 


fdh/dx \ 

\dip/dx J 


= dh/dx, (25) 


where we used Us > 0 and Ea > 0. Eq. (25) is a contradic¬ 
tion to Eq. (23). Hence, Ri < oo and thus et > 0 (water is 
entrained if flow turbulence is present). 

3.1.5. Showing < oo 

R^ < 00 can be shown by presuming R^ = oo and ar¬ 
riving at a contradiction. Inserting this presumption in the 
limit a: —>■ 00 of Eqs. (3) and (4) using Eq. (1) yields. 


dh/dx = sgn(l — Ri] c 


|^ = _3gn(l-R4c 


(26) 

(27) 


R^ = dh/dx, (28) 



Figure 1. Comparison between the downstream evolu¬ 
tions of R^ and Ri computed with the steady TEM using 
the parameter values and ignition values specified in F85 
(see Table 2) and their analytically derived asymptotic 
profiles (Eqs. (31) and (32)). 



Figure 2. Comparison between the downstream evolu¬ 
tions of U and tp computed with the steady TEM using 
the parameter values and ignition values specified in F85 
(see Table 2) and their analytically derived asymptotic 
profiles (Eqs. (35) and (33)). 
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If R-4, — 0, we cannot separate /fi) = R^ "ip/h, as we 


haviors, 


did in Eq. (24). However, in this case, we can calculate 


0 


g / dh/dx 

\EsVs-roVsC J’ 


(29) 


V’oo 

hoo 


EsVsX, 

R^x = 0.75e^x, 



(33) 

(34) 

gEsVs{8S - 5ew) \ i 

10E2 + 8Cd J ® ’ 


(35) 


from Eqs. (14) and (16), where we used > 0 and Eg > 0. 


where the subscript ’oo’ indicates the asymptotic behavior. 
We note that Eq. (33) is consistent with Eq. (16) since 


The only way in w hich t he right hand side of Eq. (29) 
vanishes is through dh/dx = 0 since 0 < < oo. Hence 

Eq. (28) is also fulfilled if R^ = 0. 

Now we calculate Ri in a similar manner and rearrange 
for using 0 < Ri < oo. This yields 


Ri = gtp/U'^ = g 


dtp/dx 


W^dU/dx 


^dU/dx = 


1 , 

3® 



R'fjj 

^dU/dx 



(30) 


—t ^ \ 1/3 

tp^ AEsVs / 10eZ + 8CD \ -i 

Coo = , „ = -— a: 3 ^ (36) 

hooUoo 3e™ \gEaVs(8S — 5eZ) J 

which vanishes in the limit a; —>■ oo. We further note that the 
derived asymptotic profiles (Eqs. (31-36)) are in agreement 
with numerical steady TEM simulations of self-accelerating 
TCs, as can be seen in Figs. 1-3, which supports the cor¬ 
rectness of our derivations. In these figures, we compare the 
downstream evolutions of R^, Ri, tp, h, U, and C computed 
with the steady TEM using the parameter values and igni¬ 
tion values specified in F85 (see Table 2 in Section 4) and 
their analytically derived asymptotic profiles (Eqs. (31-36)). 


where we used 0 < g < oo and Eq. (13). From u dx ~ 3.2. Existence criteria for self-accelerating TCs 


Cw (see Eqs. (3) and (4)) and Eqs. (1), (4), (28), and (30) 
then further follows 


Rip 

Ri 


o.Lssz:, 
lOel^ -j- 8Cd 
8S — 5eZ 


(31) 

(32) 


where we used 0 < Ri < oo, which implies S' > |el^, and 
thus self-accelerating TCs do not exist if S < Finally, 


from Eqs. (8), (16), (28), (30), (31), and (32), one obtains 


that U, tp, and h must follow the following asymptotic be- 



Figure 3. Comparison between the downstream evolu¬ 
tions of h and C computed with the steady TEM using 
the parameter values and ignition values specified in F85 
(see Table 2) and their analytically derived asymptotic 
profiles (Eqs. (34) and (36)). 


In this part of the proof, we show that < 1 is a 
necessary and sufficient condition for the existence of self- 
accelerating TCs using the results of Section 3.1. Since the 
densimetric Froude number for gravity currents is defined 
as Fr = \l\J~Ri [Kostic and Parker, 2006], Ri < 1 corre¬ 
sponds to supercritical flow. Our strategy to show Ri < 1 is 
as follows. First, we slightly modify the definition of Rip in 
Eq. (5) in a manner that ensures that the asymptotic self- 
accelerating solution of Eqs. (3-5), given by {hoo,Uoo,tpoo) 
(see Eqs. (33-35), becomes an exact solution of the mod¬ 
ified problem. Since the original problem becomes arbi¬ 
trarily close to the modified problem for self-accelerating 
TCs sufficiently far downstream, it can be considered as a 
small perturbation of the modified problem. Hence, a self- 
accelerating solution of the original problem exists if and 
only if {hao,Uoo,tpoo) is stable against small perturbations 
of the modified problem, which we show to be equivalent to 
Ri < 1. 

3.2.1. The modified problem 

We define the modified value of Rip, indicated by a tilde, 
as 

Rip = ^sVah/tp. (37) 

This definition ensures that Rip becomes arbitrarily close 
to Rip (see Eq. (13)) for self-accelerating TCs because Ea 
becomes arbitrarily close to 0 < Es < oo, while Varo/U 
becomes arbitrarily small due to l} = oo, Va < oo, and 
rt < oo. Using Eqs. (1) and (37), Eqs. (3-5) are redefined 


as 




dh 

— RiS -1- Bn, 

.i2-0.5Ri) +Cd + ^RiRip 

(38) 

dx 


1-Ri 

dU 

U RiS — Cixil + 0.5Ri) — Cd — ^RiRip 

(39) 

dx 

h 

1- Ri 

dtp 

dx 



(40) 

3.2.2. 

Stability 




It can be easily verified that {hoo,Uoo,tpao), given by 
Eqs. (33-35), is an exact solution of Eqs. (38-40). 
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We now determine the eigenvalues (A) of the Jacobi 
matrix {d{dh/dx,dU/dx,d'tp/dx)/d{h,U,ip)) evaluated at 
{hoojUtxjipoo), reading 


d{dh/dx,dU/dx,d'tp/dx) 

dih,u,^) 


(hoc , Ucc 5 ^ OO ) '— 0, 


(41) 


where I is the identity matrix, and | • | denotes the de¬ 
terminant. If the real parts of all (some) eigenvalues 
are negative (positive), small perturbations from the so¬ 
lution {hoo,Uoo,tpoo) decline (grow) with x, which means 
{hoo, Uoo, 4’<x, ) is stable (unstable). However, if the real parts 
of some of these eigenvalues vanish, one has to take a closer 
look. Indeed, one of the eigenvalues we obtain from Eq. (41) 
vanishes (Ai = 0) because d{dtp/dx)/d{h,U,'tjj) = 0 (see 
Eq. (40)). However, this does not make the solution un¬ 
stable because d{d^p/dx)/d{h,U,'il!) — 0 ensures that any 
deviation from tpoa remains constant with x. The other two 
eigenvalues (A 2 / 3 ) we obtain from Eq. (41) read 


^ 2/3 — 


-A ± VA2 - B 


(42) 


16/i(l — Ri) 

A = 24Cd + e^(28 - 

B = 24el^(l — lti){8CD — 8e^Ri{2 + Ri) + el^(10 -I- iti)), 


where we used e'^ = de^/dRi < 0, the definitions of Ri and 
R^, and Eqs. (31-35). One can see that, if Ri < 1, the 
real parts of both A 2 and A 3 are negati ve^ and the solution 
{hac,Uoo,'4>oo) thus stable. However, if Ri > 1, the eigen¬ 
value A 3 = {—A—^/A^ — B)/{16h{l —Ri)) is always positive 
due to H < 0 and the solution (hoo, Uoo, V’oo) thus unsta¬ 
ble. Furthermore, if Ri = 1, the solution is still unstable 
because small perturbations for which Ri approaches unity 
from above still result in a positive eigenvalue A 3 . Hence, 

iti < 1 (43) 


is a criterion for the stability and thus existence of self- 
accelerating TCs simulated with the steady TEM. This cri¬ 
terion leads to another criterion for S, reading 

1 C IK 

S>—Ei, + CD>—e^\m=i + CD, (44) 

O O 


where we used de^/dRi < 0 and Eq. (32). 


3.3. Sufficient conditions for Cn or alternatively S 

In this part of the proof, we presume < 0 and will 
arrive at necessary conditions for Cd and S using the results 
of Sections 3.1 and 3.2. Hence, the opposite conditions for 
Cd or alt erna tively S are sufficient for AR > 0. First, we 
show that Ai^ < 0 is never fulfilled since it leads t o a co^ tra- 
diction: In fact, AJ? < 0 necessarily implies that dk/dx < 0 
(see Eq. 6) and thus k = —00 (since k = fee, where kc is 
an arbitrary finite value, would imply dk/dx = 0), which is 
a contradiction since the Heaviside function 0(fc) in Eq. (6) 
ensures that AE = dk/dx = 0 if fc = 0. Hence, we only 
have^ to consider the case A^ = 0. In this case, we obtain 
k/U^ = 0 because, if k — 00 , 




= 0 , 


(45) 


where we used e™ > 0, Eqs. (6), (30), and (31) and 

rHospital’s rule [Chatterjee, 2012], and if fc < 00 , k/U^ = 0 
anyways. Now we perform the limit a: —> 00 on Eq. (6) us- 
ing 1 } = 00 , E^. (1), Ai^ = 0, and RiR^ — RiR^f, due to 
Ri < 00 and R^ < 00 , yielding 


= A^ = : 


0, -lti) + CD- 


(46) 


where we further used eZ < 00 , Ri < 00 , and Vs < 00 g nd 
thus RiVa/ U ~ 0, and that ke-w/U^ = 0 and /U^ = 

/3(fc/[/^)®^^ = 0 due to < 00 , < 00 , and Eq. (45). 

The maximum function in Eq. (46) occurs because, if 
0.5eZ{l — Ri) + Cd — 0.5Ri R^ is negative, k will continue 
to decrease until it vanishes, in which case 0(fc) and thus 
AE calculated by Eq. (6) also vanish. Eq. (46) provides a 
lower limit for ^RiR^, reading 

0.5e;):(l-)^)-bC'D < 0 . 5 )^)^. (47) 

Now we insert Eq. (31) into Eq. (47) and rearrange for Ri, 
yielding 


Ri > 


4 

7 




(48) 


where we used Cd > 0 and 0 < e™ < 00 . 

3.3.1. A sufficient condition for Cu 

Rearranging Eq. (48) for Cd and using de-w/dRi < 0 and 
Eq. (43) yields 


Cd ^ ^ — — C*_Dmin- ('^^) 

00 7 

Eq. (49) is a condition which must be fulfilled if AJ^ < 0. 
This means, if Cd > Gomin, self-acc eler ating TCs simulated 
with the steady TEM always fulfill Ai^ > 0. 

Finally, rearranging Eq. (32) for S and using Eq. (48) 
yields 


S — TrCiij ~\~ ■ 


+ Cd 


+ Cd 


45 

< -Cu 

- 16 


Ri 
= 5'm 


- 8®'“+4®“'ir+2Gr3 


where we used de^/dRi < 0 and that 


d 


+ Cd 


dCD \EZ + 2 Cd J 
and thus 


U- 


+ Cd 


6w -t- 2Cd 


< 


+ 0 


:-b0 


{EZ + 2Cd)^ 


< 0 


(50) 


(51) 


(52) 


Eq. (50) is a condition which must be fulfilled if AJ^ < 0. 
This means, if S' > Smin, self-acce lera ting TCs simulated 
with the steady TEM always fulfill Ai^ > 0. 


3.4. Are Cd > Comin or S > Smin? 

Now we argue that at least one of the conditions Cd > 
Gomin and S > Smin is virtually always fulfilled. Ac¬ 
cording to Eqs. (49) and (50), the values of Gomin and 
Smin depends on the empirical relationship used to calcu¬ 
late Cu,. Standard relationships for e™ are, for instance, 
Eq. (9), for which Gomin = 0.00097 and Smin = 0.0073, and 
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Table 2. Parameter values and ignition values specified in 
F85 and P86. The values in brackets correspond to the mod¬ 
ified ignition values which we used to obtain self-accelerating 
TCs for the parameter values specified in P86. 



Fukushima et al. [1985J 

Parker et al. [1986J 

s 

0.08 

0.05 

Da [mm] 

0.15 

0.1 

Ds [m/s] 

0.0165 

0.0084 

y [mVs] 

10-® 

10-® 

Pa [kg/m3] 

2650 

2650 

pn, [kg/m®] 

1000 

1000 

Cd 

0.004 

0.004 

h(0) [m] 

3 

2 (1) 

U(0) [m/s] 

1.24 

0.652 (0.699) 

m [mVs] 

0.019 

0.0038 (0.0047) 


e™ = 0.075/Vl + 718-Ri2-^ [Parker et al., 1987], for which 
Comin = 0.00205 and 5'min = 0.0154. However, recent 
state of the art simulations of suspended sediment trans¬ 
port [Schmeeckle, 2014] (see Section 5.2) as well as the ma¬ 
jority of experimental data [Bradford and Katopodes, 1999] 
indicate that realistic values of Cd are significantly larger. 
Hence , in realistic simulations self-accelerating TCs always 
fulfill A7^ > 0. Even if an unrealistically small value for Cd 
is used, one would still have to consider bed slopes which 
are much smaller than those one typically uses for TCs. 
For instance, F85 used Cd = 0.004 and S = 0.08. Each 
of these values alone is actually sufficient to ensure that 
thei r sim ulated self-accelerating TCs were ph ysic ally realis¬ 
tic (Ai^ > 0). The fact that F85 reported Ai^ < 0 thus 
means that there was most likely an error in the numeri¬ 
cal computations. Indeed, we show in the following that, 
depending on the ignition values, numerical simulations us¬ 
ing the same parameter values and empirical relations as 
reported by F85 and P86 result either in self-accelerating 
TCs which fulfill Ai^ > 0 or in TCs which fulfill ij < oo 
and thus are not self-accelerating. 

4. Numerical simulations using the steady 
TEM 

In this section, we first report simulations of self- 
accelerating TCs using the steady TEM and exactly the 



Figure 4. TCs simulated with the steady TEM for the 
parameter values and ignition values specified in F85 and 
P86 (see Table 2). The solid lines show the downstream 
profiles of U and the dashed lines the downstream pro¬ 
files of i/j. The star labels the case in which the ignition 
values used in P86 are modified to those given in brackets 
in Table 2. 


same empirical relations (Eqs. (9-11)), physical parame¬ 
ters and ignition values (see Table 2) as F85 and P86 in 
Section 4.1. 

Afterwards in Section 4.2, we present a new method to 
obtain ignition values, which even results in self-accelerating 
TCs if Ri is very close to unity. 

In all our simulations, the model equations were inte¬ 
grated using the Runge-Kutta method, and we confirmed 
that reducing the spatial integration step did not signifi¬ 
cantly change the results. 

4.1. Repetition of the original simulations by F85 
and P86 

Fig. 4 shows the downstream evolutions of U (solid lines) 
and tp (dashed lines) for the TCs simulated with the steady 
TEM. As can be seen, while the parameter values and igni¬ 
tion values specified in F85 result in self-accelerating TCs, 
those specified in P86 do not since U and ip are decreasing 
downstream for x > 212m. In order to avoid confusion, it 
is very important to realize here that P86 called this TC 
“self-accelerating” because of its accelerating behavior for 
X < 212m, while it is not self-accelerating according to our 
definition. Moreover, these authors found that the decrease 
of U and ip for x > 212m coincides with AE < 0. How¬ 
ever, it is not surprising at all that decelerating TCs loose 
turbulent kinetic energy when moving downstream. In fact, 
also the steady FEM produces decelerating TCs if the initial 
conditions are chosen appropriately [Parker et al., 1986], and 
these currents should generally also loose turbulent kinetic 
energy when moving downstream. There is no qualitative 
difference between the steady TEM and FEM in this re¬ 
gard. Hence, in order to disproof the claim of P86 tha t the 
steady TEM cannot produce physically realistic (AJ^ > 0) 
self-accelerating TCs, we only have to show that, for the 
same physical parameters {S, Da, Va, u), there are igni¬ 
tion values which res ult^ in self-accelerating TCs since these 
automatically fulfill A7^ > 0 (see our proof in Section 3). 
Indeed, by changing the ignition values (see the values in 
brackets in Table 2), the parameter values specified in P86 
result in self-accelerating TCs (the case labeled by a star in 
the legend of Fig. 4). To be consistent, we used the same 
procedure as P86 (the method by Parker ]1982]) to obtain 
the modified ignition values, which is defined in the follow¬ 
ing. First, one chooses arbitrarily the upstream boundary 



Figure 5. TCs simulated with the steady TEM for the 
parameter values and ignition values specified in F85 and 
for the parameter values specified in P86 with the modi¬ 
fied ignition values (see Table 2), indicated by a star. The 
cases p = 1.5 and p = 0.5 refer to the use of a modified 
bed sediment erosion rate relation (Eq. (55)). 
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condition h{0). Then one defines a quasi-equilibrium state 
by {dU/dx){0) = {dipjdx){Q) — 0, but {dh/dx){0) / 0. Us¬ 
ing this Eqs. (4) and (5) can be solved for U{0) and ip{0). 
While P86 originally chose h{0) = 2m, from which they 
computed U{0) = 0.652m/s and ip{0) = 0.0038m^/s, the 
modified value reads h{0) = Im, from which one computes 
U{0) = 0.699m/s and ip{0) = 0.0047m^/s (see Table 2). 

Note that the fact that the original ignition values spec¬ 
ified by P86 result in a decelerating TC shows that the 
method by Parker [1982] to obtain the ignition values does 
not always work very well. The reason is that this method 
does not necessarily ensure that the ignition state is suffi¬ 
ciently close to the asymptotic solution {hoo,ipoo,Uoo)- In 
Section 4.2 we will provide an improved method to obtain 
the ignition values that even works in cases in which t he 
method by Parker [1982] does not at all, namely when Ri 
is close to unity. 

It remains to show that the self-accelerating TCs shown 
in Fig. 4 (i.e., the P85 current and the modified P86 cur¬ 
rent), but not the original P86 current (which is not self- 
accelerating, but decelerating), fulfill > 0. In order to 
do so, we note that, if 


ASmod = U’’ ( -e^(l — Ri) + Cd — ~ —RiR^ 


(53) 


is smaller than zero, then, due to Eqs. (1) and (6) and 
e™ > 0, /? > 0, and k > 0, 


AE - ^3 ASmod ^2 


/3fc3/2 

t/3 


(54) 


is also smaller than zero. It follows that AEmod > 0 is a 
condition which must be fulfilled ever after a hnite distance 
downstream in orde r for self-accelerating TCs to be phys¬ 
ically realistic (Kii >0). In fact, this criterion was used 
by F85 and P86 to determine whether self-accelerating TCs 
are physically realistic. Indeed, in accordance with our an¬ 
alytical proof in Section 3, Fig. 5 shows that AEmod > 0 
is fulfilled. It shows the downstream profiles of AiJniod for 
cases in Table 2 that resulted in self-accelerating TCs. It can 
be seen that the qualitative behavior of Aifniod is exactly 
opposite to the one described in F85. Instead of AEmod be¬ 
ing positive in the beginning and more and more negative 
later on downstream, as described by F85, AUmod is nega¬ 
tive in the beginning and more and more positive later on 
downstream. This strongly suggests an error in the early 
computations by F85. 

Fig. 5 also shows the downstream profiles of ASmod for 
the parameter values and ignition values specified in F85, 
but with a modified bed sediment erosion rate relation. 


Es 


0.3p ^ > 13.2 'I 

3p X 10-i2^i°(l - 5/^) 5 < ^ < 13.2 I, (55) 

0 Z <h \ 


where p > 0 is a constant factor (the case p = 1 is identical 
to Eq. (11)). It can be seen that the larger is the value of p 
and thus Es the larger is AEmod- This is again in contrast 
to the claim of F85 and P86 that their ostensible failure of 
the steady TEM is due to strong erosion of bed sediment. 
According to their argument, stronger erosion should lead to 
more energy spent in eroding and suspending bed sediment 
and thus to smaller values of AEmod. However, one must 
also take into account that the eroded bed sediment has a 
potential energy which is converted into turbulent kinetic 
energy when moving downslope. This conversion would not 
take place if the sediment continued to rest at the top of the 
sediment bed. It can be seen in Fig. 5 that this manner of 


additional production of turbulent kinetic energy more than 
compensates the additional energy loss due to the stronger 
erosion of bed sediment. Indeed, as we show in Section 5.1, 
AUmod is proportional to Es- 


4.2. Ignition of self-accelerating TCs 

In S^tion 3.2, we showed that self-accelerating TCs ex¬ 
ist for Ri < 1 when simulated with the steady TEM. This 
might seem quite surprising because the condition Ri > 0.25 
is generally used to estimate whether the turbulent mixing 
becomes insufficient to overcome density layering in physi¬ 
cal environments [Prandle, 2009]. For this reason, we show 
here that numerical simulations with the steady TEM, in¬ 
deed, result in self-acceleratin g T Cs if the ignition values are 
chosen appropriately, even if Ri is very close to unity. 

In order for {h, U, ip) to converge against the asymptotic 
solution (hoc, Uoa,ipao), the upstream boundary values of Ri 
and R^ should not deviate too strongly from their asymp¬ 
totically constant values Ri and R^ since we only showed 
stability against sufficiently small deviations in Section 3.2. 
We thus propose to use the following ignition values 


Ri{0) = Ri 


R^ ( 0 ) = R^ 


5a Cd b 

1 ^ ^ ~ 2 _ 

// 5a Cd b\^ 5a CDb 

~ 2 ) 

3a 

4(6 + Ri(0))’ 


(56) 

(57) 


where a = 0.00153, p = 0.02fl4, and the right hand sides are 
the solutions for R^ and Ri of Eqs. (31) and (32), respec¬ 
tively, when Cu, is calculated by Eq. (9). For a given value 
of h{0), Eqs. (56) and (57) can be solved to obtain the up¬ 
stream boundary values t/(0) and ip{0) using Ri = gip/U^ 
and Eq. (5). Thereby larger values of h{0) correspond to 
larger values of U{0) and ip{0) and are thus relatively closer 
to the asymptotic solution. Hence, a sufficiently large value 
of h{0) ensures that the ignition values calculated in the 
manner above always leads to self-accelerating TCs. 

Fig. 6 shows the downstream evolutions of U, ip, and h 
computed with the steady TEM using the parameter val¬ 
ues as specffied in P86 (see Table 2), but with S = 0.0069 



Figure 6. Downstream evolutions of 17, ip, and h (inset: 
Ri) computed with the steady TEM using the parameter 
values specified in P86 (see Table 2), but with S = 0.0069 
instead of S' = 0.05. The ignition values are h{0) = 12m, 

Ri{0) = iti = 0.99, and R,i,{Q) = R^ = 0.00114 (corre¬ 
sponding to U{0) = 5.58m/s and ip{0) = 10.64m^/s). 
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instead of S' = 0.05. This modified value of S corresponds 
to Ri = 0.99 (from Eq. (32)). The upstream depth of 
the TC is set to h(0) = 12m. Then 17(0) = 5.58m/s and 
'tp{0) = 10.64m^/s are obtained from the conditions Ri{0) = 
Ri = 0.99 (from Eq. (56)) and R^{0) = R^ = 0.00114 (from 
Eq. (57)). It can be seen that our improved definition of 
the ignition values even leads to self-accelerating TCs when 
Ri is very close to unity, as analytically predicted. 

We find that larger values of h{0) also result in self- 
accelerating TCs, while significantly smaller values of h{0) 
do not, in agreement with our analytical prediction that 
self-accelerating TCs always exist if h{0) is larger than a 
certain minimal value. Indeed, for all values h{0), Ri ini¬ 
tially increase s d ownstream before it decreases and con¬ 
verges against Ri (see inset of Eig. 6). If Ri exceeds unity in 
its initial increase, the TC rapidly dies, and this always hap¬ 
pened if a significantly smaller value than 12m had been cho¬ 
sen for h{0). The upstream boundary condition h{0) = 12m 
is, however, just sufficient to ensure that Ri does not ex¬ 
ceed unity in its initial increase (see inset of Eig. 6). We 
further confirmed our analytical jrrediction in Section 3.2 
that any value of S leading to Ri > 1 never resulted in 
self-accelerating TCs. 

We would like to emphasize that the method described 
above to determine the ignition values is a major improve¬ 
ment over the method given by P86, in which h{0) is 
fixed, and (7(0) and ip{0) are obtained from {dU/dx){0) = 
{d'(p/dx){G) = 0. If fact, for the parameter values specified 
above, which correspond to Ri = 0.99, it is impossible to ob¬ 
tain self-accelerating TCs with the method by P86 regardless 
of the value we chose for h{0). Even for parameter values 
for which Ri is significantly smaller unity, it is remains un¬ 
certain for which values of /i(0) this method yields ignition. 
Eor instance, for the parameter values as specified in P86 
(see Table 2), h{0) = Im yields ignition and h{0) — 2m does 
not. 

Another major improvement of the method described 
above to determine the ignition values is that h, U, and 
tp follow the computed asymptotic behaviors (Eqs. (33-35)) 
already at the upstream boundary (see Eig. 6), while they 
otherwise might require very long distances to do so and thus 
be extraordinarily large. The fact that h, (7, ip are within 
a realistic range at the upstream boundary when it is de¬ 
termined using this improved method (e.g., for S — 0.05, 
a corresponding ignition condition would be h{0) — 0.5m, 
U = 0.87m/s, and ip = O.Olm^/s) shows that the asymptotic 
analysis we performed in Section 3 is not just a mathemat¬ 
ical exercise, but has real-world relevance. 


5. Discussion 

This section contains two parts: a mathematical expla¬ 
nation for the increase of ARmod with Eg in Eig. 5 in Sec¬ 
tion 5.1 and a discussion of the question whether the TEM 
or EEM is more realistic in Section 5.2. 


5.1. Relation between bed erosion and net tnrbulent 
kinetic energy prodnction 

In this section, we discuss the consequences of Ai^ > 0 for 
the turbulent kinetic energy k. This leads to a relation be¬ 
tween the bed sediment erosion rate and the net production 
rate of turbulent kinetic energy which explains the results 
shown in Eig. 5. In order to do so, we first show that 
xdk/dx = 0 if k = fcc, where kc is an arbitrary finite value. 
In fact, this follows from 



-r^dk, dk 

k + X— = fcc + X — , 

dx dx 


(58) 


where we used I’Hospital’s rule [Chatterjee, 2012]. Hence, 
also AE = {h/U^)dk/dx ~ x^^^dkjdx ^~*°°> 0 ( from 
Eqs. (34) and (35)) if = k^, which c ontrad icts Aik > 0. 

This means that fc = oo and thus k/U^ = 2Akje^, 
which follows from Eq. (45) and 0 < ejj < oo. Hence, 

e^k/U^ = 2Ak and /U^ = 2~$(Kk . With 
this knowledge, one can now calculate the limit a: —>■ oo of 
Eq. (6), analogous to what we did in Eq. (46) for the case 
AR = 0, yielding 

iAk + 2~$(Ak/etf/^ = ^Et{l-lk) + CD - 

(59) 

Eq. (59) can be solved for Ak and implies that 'Kk < oo 
since 0 < ejj < oo, Ri < oo, Cd < oo, and R^< oo. 
Hence, once one has determined the value of from 
Eq. (59), one can compute the asymptotic profile of k from 
dk/dx = {U^/h)AE and Eqs. (34) and (35), reading 



R^ 



(60) 


when inserting the values of R^ and Ri (Eqs. (31) and (32)). 

We wish to emphasize that the asymptotic profiles of h, 
Ip, U, and k of self-accelerating TCs computed with the 
steady TEM (Eqs. (31-35) and (60)) are identical to the 
same p rofiles computed by the FEM. if Cd is replaced by 

ak/U^ and one assumes 0 < ak/U^ < oo. This is be¬ 
cause the derivations in Section 3.1 remain the same in 
this case (see our statement in the first paragraph of Sec¬ 
tion 3.1) and ak[ U^ = 2aAk/e^ (see Eq. (45)), from which 
follows 0 < AR < oo. Interestingly, the physical mean¬ 
ing of the quantity a.k/U'^ in the FEM is exactly that of 
the bed drag coefficient [Parker et ai, 1986]. This means 
that physically relevant self-accelerating TCs computed with 
the FEM ( those with finite, positive bed drag coefficient, 

0 < ak/U^ < oo) are qualitatively identical to those com¬ 
puted with the steady TEM since only the prefactors in the 
asymptotic profiles are differen t. We note that it can prob¬ 
ably be shown that 0 < ak/U^ < oo for all self-accelerating 
TCs computed with the FEM. 

Using the results above, we now take a look at the di¬ 
mensional net production rate of turbulent kinetic energy, 
given by Rnet = AE + Cui Uk [P arker et al, 1986]. In both 
the steady TEM and FEM, Ak i s indep endent of the ero¬ 
sion rate Eg since et, , Cd (or ak/U in the FEM), R^p, 
and Ri are independent of Eg. This mean that the asymp¬ 
totic dependency of Rnet on Eg is entirely incorporated in 
and Uk. From Eqs. (35) and (60), we thus learn that all 
contributions to production and dissipation of turbulent ki- 
ne tic energy and thus Rnet are asymptotically proportional 
to Eg. In fact, not only the dissipation due to erosion of bed 
sediment (0.5RiR,/,) is asymptotically amplified by Eg, as 
argued by F85 and P86, but also the dissipation due to wa¬ 
ter entrainment (0.5eu,Ri), viscous dissipation {/3k^^^), and 
turbulent kinetic energy production {0.5emU^ -|-m*R). This 
eventually explains the increase of ARmod with Eg in Fig. 5, 
which since ARmod = Rnet + and thus ARmod is also 

asymptotically proportional to Eg. 


5.2. An attempt to compare the TEM with the FEM 

In this section, we attempt to compare the physical re¬ 
alism of the TEM with that of the FEM. While the pre¬ 
vious parts of the paper dealt with the steady TEM and 
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Figure 7. Comparison between the TEM and FEM clo¬ 
sures for the bed shear stress. The blue circles correspond 
to 55Cd = 55u^/U^ (see Eq. (l))and the red triangles 
to a = (see Eq. (2)), where w», U, and k are ob¬ 
tained from the suspended load data plotted in Figs. 4 
and 6-8 in Schmeeckle [2014]. Constant values of 55Cd 
and a correspond to perfect agreement with the respec¬ 
tive closures. Cd has been multiplied by factor 55 to 
make visual comparison between the closures easier. 


FEM versions proposed by F85 and P86, we now attempt 
to make a more general assessment since a large number of 
improved TEM and FEM versions have been proposed since 
these original studies were published. Because of this, the 
most meaningful way to compare the TEM with the FEM, in 
our opinion, is to evaluate the equation which all TEM and 
FEM versions, respectively, have in common and in which 
the TEM differs from the FEM: the closure for the bed shear 
stress (Eqs. (1) and (2), respectively). In order to do so, we 
first briefly reiterate the assumptions behind these closures. 

On the one hand, Eq. (1) follows from the idea that the 
fluid shear stress at the bed (p^ul) describes the streamwise 
component of the force applied by the fluid on the station¬ 
ary bed per unit area. The main streamwise fluid force is 
the mean drag force, which is proportional to the square of 
the mean local flow velocity (u). Hence, this assumption 
yields u* oc u(zi,)^, where zi, denotes the vertical location 
of the top of the sediment bed. u(zi,) is then assumed to 
be roughly proportional to U = j u{z)dz (the height- 

averaged flow velocity), which eventually yields uj oc U^. 
On the other hand, Eq. (2) follows from the definition of the 
Reynolds stress (u« = —u'v'{zb))- —u'v'{zb) is first assumed 
to be proportional to 0.5u'^(z6) because both are turbulent 
correlations of dimension velocity square. Then 0.5u'2(z6) 
is assumed to be proportional to k — 0.5-^ u'^{z)dz 

(i.e., the height-averaged value of 0.5u'2), yielding uj oc k. 

In the following, we compare both closures with the sim¬ 
ulations of turbulent capacity sediment transport (using 
a coupled Large Eddy and Discrete Element Model) by 
Schmeeckle [2014]. These simulations belong to the most 
realistic turbulent sediment transport simulations in the lit¬ 
erature because they consider several layers of the particle 
bed at the scale of the particle, including particle-particle in¬ 
teractions as well as momentum extraction from flow due to 
drag on the particles. Schmeeckle [2014] simulated bedload 
and suspended load. However, only the suspended load sim¬ 
ulations are of interest for us since both the TEM and FEM 
assume that the instantaneous and thus average horizon¬ 
tal particle and flow velocities are the same and that there 


is thus no horizontal fluid drag term in the horizontal mo¬ 
mentum balance of the fluid. In other words, the fluid shear 
stress, whose vertical gradient appears in this horizontal mo¬ 
mentum balance, is assumed to be undisturbed by the pres¬ 
ence of transported particles. However, as shown in Fig. 10 
in Schmeeckle [2014], this assumption is only fulfilled in the 
upper parts of the flow (zjh > 0.3, where h in Schmeeckle 
[2014] is the simulation height). Hence, the “bed” height 
{zb) in the TEM and FEM is actually the height above which 
the local fluid shear stress is undisturbed by the presence of 
transported particles, which implies Zb ~ 0.3/i for the sus¬ 
pended load simulations by Schmeeckle [2014]. Note that 
u*(l — 0.3) in Schmeeckle [2014] thus corresponds to m* in 
this manuscript. Using this value of Zb, we used the sus¬ 
pended load data plotted in Figs. 4 and 6-8 in Schmeeckle 
[2014] to compute Co = and a = uljk (see Eqs. (1) 

and (2)). The values of 55(70 and a obtained in this way 
are plotted in our Fig. 7. 

It is important to note that we excluded the simulation 
with strongest suspended transport from this figure because, 
exclusively in this simulation, there is significant flow speed 
in the entire simulation domain (see Fig. 4 in Schmeeckle 
[2014]). This implies a reduction of the flow resistance and 
thus Cd which can be attributed to the finite size of the 
simulated system. To explain this, let us imagine, we ex¬ 
tend the simulation domain by adding a layer of particles 
below 2 = 0 and moving the lower simulation wall to the 
bottom of this added particle layer. Because the flow speed 
near the added particle layer is significant in the simulation 
with strongest suspended transport, such an extension of the 
simulation domain would lead to an increase of the overall 
horizontal drag on the particles and thus to increasing flow 
resistance, while such an extension of the simulation domain 
would have nearly no effect in the other simulations due to 
zero flow speed near 2 = 0. 

It can be seen in Fig. 7 that the simulations by 
Schmeeckle [2014] seem to slightly support the use of the 
TEM over the FEM closure since Cd varies slightly less with 
M* than a, which seems to slightly increase with w*. We be¬ 
lieve that this increase is the result of turbulence damping 
due to density stratification, as we explain in the follow¬ 
ing. On the one hand, in order to keep the stratification 
stable, the flow must exert vertical drag forces on the par¬ 
ticles which on average exactly compensate the submerged 
gravity forces. However, through these vertical drag forces, 
the flow loses turbulent kinetic energy (fc). With increasing 
suspended load, the concentration and thus the submerged 
weight of the particles increases, resulting in a decrease of k. 
On the other hand, the vertical fluid shear stress profile and 
thus u* in the TEM and FEM are by definition undisturbed 
by the presence of transported particles (since the TEM and 
FEM assume that the instantaneous and thus average hori¬ 
zontal particle and flow velocities are the same) and thus not 
influenced by density stratification. Hence, ut/k increases 
with u*. 

We wish to emphasize, since the increase of a with u* is 
quite small, our comparison is just a first clue in favor of 
the TEM, which needs to be further supported by data in 
the future. Also, it is important to check in the future how 
both Cd and a depend on flow parameters which remained 
constant in the simulations by Schmeeckle [2014], such as 
the particle Reynolds number {Rcp). 

6. Conclusions 

This study re-examines the steady three-equation model 
(TEM) for turbidity currents (TCs) by Fukushima et al. 
[1985] (F85) and Parker et al. [1986] (P86) by analytical and 
numerical means and compares the TEM and four-equation 
model (FEM) closures with predictions of recent numerical 
simulations. The following conclusions can be drawn from 
this study: 
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1. Self-accelerating TCs simulated with the steady TEM 
by (F85) and (P86) never violate the turbulent kinetic en¬ 
ergy balance if a realistic value for the bed drag coefficient 
(Cd) is used (see Sections 3 and 4), which is in contrast 
to the nearly three decades old scientific consensus on that 
matter. 

2. The asymptotic behaviors of self-accelerating turbid¬ 
ity currents have been analytically calculated in Sections 3.1 
and 5.1 (see Eqs. (34-36) and (60)). 

3. It is not necessary to limit the bed erosion rate (Es) 
to allow for self-accelerating TCs, which was the motiva¬ 
tion for the FEM, since the net production rate of turbu¬ 
lent kinetic energy is asymptotically proportional to Es (see 
Section 5.1), even though turbulent kinetic energy is spend 
when suspending bed material. The physical reason behind 
this counter-intuitive behavior is that eroded bed sediment 
increases the sediment concentration and thus potential en¬ 
ergy of the TC, which is then converted into turbulent ki¬ 
netic energy downslope. 

4. The steady TEM investigated in this paper has nu¬ 
merically stable self-accelerating solutions if and only if the 
Richardson number (Ri) is smaller than unity (supercriti¬ 
cal flow) in the asymptotic limit (Ri = lima;_nx) Ri < I, see 
Section 3.2), which can be calculated by Eq. (32). This 
condition is equivalent to the condition that the bed slope 
{S) must be larger than a critical value (see Eq. (44)). 

5. A novel method to determine the ignition values has 
been proposed (see Section 4.2). This method always leads 
to self-accelerating TCs if Ri < 1 in simulations with the 
steady TEM investigated in this paper, which is a major 
improvement over the method by Parker [1982], which of¬ 
ten fails to do so. 

6 . The TEM and EEM closures for the bed shear stress 
are compared with state of the art simulations of suspended 
sediment transport using a coupled Large Eddy and Dis¬ 
crete Element Model (see Section 5.2). These simulations 
suggests that the TEM closure (Eq. (1) performs slightly 
better that the FEM closure (Eq. (2). 

It is important to mention that most if not all of the 
conclusions summarized above can be easily generalized to 
more modern versions than the E85 and P86 version of the 
TEM investigated in this paper. Eor instance, unsteady 
self-accelerating TC solutions can be treated as fluctuations 
around the steady solution. Depending on the magnitude of 
these Jffictuations, the critical asymptotic Richardson num¬ 
ber {Ri) for the existence of stable solution then will be 
somewhat smaller than unity. In fact, just so small that 
fluctuations {Ri') of Ri never lead to Ri > 1, meaning 
Ri < 1 — maxRi'. Moreover, the conclusions summarized 
above indicate a strong need for studies comparing the TEM 
and FEM with each other in the future in order to assess 
which of these models is more realistic. Our study just pro¬ 
vides a first small clue in favor of the TEM, but more inves¬ 
tigations are needed. 
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